//----------------------------------------------------------------------------
// M A R I T I M E  R E S E A R C H  I N S T I T U T E  N E T H E R L A N D S
//----------------------------------------------------------------------------
// Copyright (C) 2011 - MARIN - All rights reserved - http://www.marin.nl
//----------------------------------------------------------------------------
// Program    : mermaid
// Module     : vbm
// File       : CudaRrbSolver.h
// Author     : M.A. de Jong
//----------------------------------------------------------------------------

#ifndef CUDARRBSOLVER_H
#define CUDARRBSOLVER_H

#include "Solver.h"
#include "Array2D.h"

#include "CUDAPACKAGE/common/globals.h"

#define CUDARRBSOLVER_BASE Solver


//----------------------------------------------------------------------------
// class CudaRrbSolver
//----------------------------------------------------------------------------
class CudaRrbSolver : public CUDARRBSOLVER_BASE
{
    typedef CUDARRBSOLVER_BASE base;

public:

    CudaRrbSolver(Array2D<REAL> *pcc,
                  Array2D<REAL> *pcs,
                  Array2D<REAL> *pcw,
                  int Nx1,
                  int Nx2,
                  REAL tolerance);

    virtual ~CudaRrbSolver();

    virtual void solve(Array2D<REAL> *pB,
                       Array2D<REAL> *pX);

    // using device pointers directly, saves memory transactions
    // this function is used when the CudaWavesComputer is enabled
    virtual void cudaSolve(REAL *bpsi,
                           REAL *psi,
                           Array2D<REAL> *pB,
                           Array2D<REAL> *pX);

private:

    /// @cond nevercall
    CudaRrbSolver(const CudaRrbSolver& m);
    CudaRrbSolver & operator =(const CudaRrbSolver& m);
    /// @endcond

    void setup();

    void cleanup();

    void preconditioner();

    void solver();

    void cg();

    void matvec(REAL *dY, REAL *dX);

    void axpy(REAL *dY, REAL *dX, REAL a, REAL b);

    void determineGrids(unsigned int Nx1, unsigned int Nx2);

    REAL dot(REAL *dY, REAL *dX);

    unsigned int m_nx[20];
    unsigned int m_ny[20];
    unsigned int m_cx[20];
    unsigned int m_cy[20];
    unsigned int m_numLevels;
    unsigned int m_remSize;

    Array2D<REAL> *m_pcc; // center stencil on host
    Array2D<REAL> *m_pcw;
    Array2D<REAL> *m_pcs;

    REAL *m_dcc; // center stencil on device
    REAL *m_dss;
    REAL *m_dww;
    REAL *m_dX;
    REAL *m_dB;
    REAL *m_dXr1r2b1b2;
    REAL *m_dYr1r2b1b2;
    REAL *m_dBr1r2b1b2;
    REAL *m_dPr1r2b1b2;
    REAL *m_dQr1r2b1b2;
    REAL *m_dRr1r2b1b2;
    REAL **m_dZr1r2b1b2;
    REAL *m_drem;
    REAL *m_hrem;

    Grid m_orig;
    Grid *m_prec;
};

inline Solver* NewCudaRrbSolver(Array2D<REAL> *pmatC, ///< central matrix
                            Array2D<REAL> *pmatS, ///< south matrix
                            Array2D<REAL> *pmatW, ///< west matrix
                            int N, ///< width
                            int M, ///< height
                            REAL tolerance) ///< solver tolerance
{
    return new CudaRrbSolver(pmatC, pmatS, pmatW, N, M, tolerance);
}

#endif // CudaRrbSolver_H
